Enfoque estadístico del aprendizaje

Trabajo práctico nº1: regresión lineal

 

Datos a utilizar:

Los datos con los que se trabajará en este TP provienen de la 3° Encuesta Mundial de Salud Escolar (EMSE) provistos por el Ministerio de Salud de la República Argentina. Esta encuesta trata sobre temas de salud y hábitos de las personas en la escuela secundaria que pueden impactar en su salud.

Los datasets que se comparten corresponden a un recorte (muestra) del dataset original, luego del tratamiento de valores atípicos e ingeniería de atributos.

Las variables incluidas son:

Variable	Descripción
record	ID de la observación
edad	Edad en años
genero	género de la persona
nivel_educativo	nivel educativo en que se encuentra la persona
altura	altura en centímetros
peso	peso en kilogramos
frecuencia_hambre_mensual	variable categórica que indica la frecuencia con la que la persona considera que pasó hambre en el último mes porque no había suficiente comida en su hogar
dias_consumo_comida_rapida	cuántos días comió en un restaurante de comida rápida en la última semana
edad_consumo_alcohol	edad en qué la persona comenzó a consumir alcohol
consumo_diario_alcohol	cantidad de tragos que la persona habitualmente toma por día
dias_actividad_fisica_semanal	cantidad de días que la persona realizó una actividad física por un total de al menos 60 minutos en la última semana
consumo_semanal_frutas	cantidad de veces que la persona consumió frutas en la última semana
consumo_semanal_verdura	cantidad de veces que la persona consumió gaseosas (al menos un vaso) en la última semana
consumo_semanal_gaseosas	cantidad de veces que la persona consumió snacks/comida salada en la última semana
consumo_semanal_snacks	cantidad de veces que la persona consumió snacks/comida salada en la última semana
consumo_semanal_comida_grasa	cantidad de veces que la persona consumió comidas altas en grasas en la última semana

Consignas:

El objetivo general del trabajo es poder crear una serie de modelos lineales para explicar y predecir el peso de los estudiantes según la información que proporciona la EMSE.

---

 Se cargan las librerías con las que se va a realizar el trabajo:

library(tidyverse)
library(tidymodels)
library(ggplot2)
library(knitr)
library(GGally)
library(robust)

options(scipen=999)

 

Importo el dataset:

encuesta_salud <- read.csv("encuesta_salud_train.csv")

1) Análisis exploratorio

Primero, observo algunos registros del dataset:

encuesta_salud %>% sample_n(5)

Observo las dimensiones del dataset:

glimpse(encuesta_salud)

Rows: 7,024
Columns: 16
$ record                        <int> 502, 26488, 31473, 14154, 36578, 53730, 11892, 35549, 14928, 11754, 40850, 11038, 32832, 31300, 32497, 6961, 39376, 30541, 18…
$ edad                          <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, 16, 14, 15, 17, 15, 14, 15, 17, 17, 16, 14, 12, 15, 17, 15, 16, 16, 16, 17, 15, 17, 1…
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "Masculino", "Masculino", "Masculino", "Femenino", "Femenino", "Masculino", "Femenino",…
$ nivel_educativo               <chr> "2do año/11vo grado nivel Polimodal o 4to año nivel Secundario", "1er año/10mo grado nivel Polimodal o 3er año nivel Secundar…
$ altura                        <int> 165, 178, 172, 170, 170, 178, 156, 163, 164, 167, 185, 146, 180, 175, 183, 165, 165, 157, 165, 170, 174, 163, 154, 164, 160, …
$ peso                          <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, 100, 33, 62, 70, 80, 60, 47, 50, 50, 70, 75, 55, 64, 59, 86, 51, 71, 49, 70, 65, 65, …
$ frecuencia_hambre_mensual     <chr> "Rara vez", "Rara vez", "Nunca", "Nunca", "Rara vez", "Nunca", "Nunca", "Nunca", "Nunca", "Nunca", "Nunca", "Rara vez", "Nunc…
$ dias_consumo_comida_rapida    <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1, 3, 0, 0, 0, 0, 1, 0, 6, 0, 1, 0, 2, 0, 2, 0, 0, 5, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0,…
$ edad_consumo_alcohol          <chr> "14 o 15 años", "7 años o menos", "Nunca tomé alcohol más que unos pocos sorbos", "14 o 15 años", "16 o 17 años", "8 o 9 años…
$ consumo_diario_alcohol        <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, 5.0, 0.0, 5.0, 0.0, 0.0, 2.0, 1.0, 0.0, 5.0, 0.0, 1.0, 1.0, 0.5, 0.0, 3.0, 3.0, 1.0, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1, 4, 0, 1, 6, 5, 7, 3, 0, 7, 5, 2, 2, 4, 2, 7, 4, 7, 0, 6, 4, 0, 3, 2, 0, 2, 2, 3, 0,…
$ consumo_semanal_frutas        <chr> "No comí frutas durante los últimos 7 días", "No comí frutas durante los últimos 7 días", "No comí frutas durante los últimos…
$ consumo_semanal_verdura       <chr> "4 a 6 veces durante los últimos 7 días", "4 a 6 veces durante los últimos 7 días", "1 vez al día", "4 o más veces al día", "…
$ consumo_semanal_gaseosas      <chr> "1 a 3 veces durante los últimos 7 días", "1 a 3 veces durante los últimos 7 días", "4 a 6 veces durante los últimos 7 días",…
$ consumo_semanal_snacks        <chr> "1 a 3 veces durante los últimos 7 días", "No comí comida salada o snacks en los últimos 7 días", "4 a 6 veces durante los úl…
$ consumo_semanal_comida_grasa  <chr> "No comí comida alta en grasa en los últimos 7 días", "4 a 6 veces durante los últimos 7 días", "No comí comida alta en grasa…

Aquí se puede ver que el dataset cuenta con 7024 observaciones para las cuales se tienen 16 variables (7 de ellas numéricas y 9 categóricas). La variable peso, que es una de las variables numéricas, es la variable a explicar con el resto.

A continuación, observamos si las variables numéricas - y fundamentalmente el peso - están correlacionadas o no. Para ello, se procede a aperturar por el género de cada una de las personas:

encuesta_salud %>%
  select(where(is.numeric), genero, -record) %>%
  ggpairs(mapping = aes(color = genero), title = "Matriz de correlaciones",
          upper = list(continuous = wrap("cor", size = 5, hjust=0.5)), legend = 25) +
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, vjust=0.5), legend.position = "bottom")

En los gráficos, se observa que existe cierta correlación entre entre el peso y la altura (0,57 para las personas de género masculino, 0,44 para las de género femenino. También hay una correlación moderada entre la edad y el peso.

encuesta_salud %>% 
  group_by(genero) %>% 
  summarise(cor = cor(altura, peso))

Luego, observamos qué categorías existen para la “frecuencia de hambre mensual”:

encuesta_salud %>%
  select(frecuencia_hambre_mensual) %>%
  table() %>%
  prop.table() %>%
  sort(decreasing = TRUE)

.
        Nunca      Rara vez Algunas veces  Casi siempre  Dato perdido       Siempre 
  0.689635535   0.206007973   0.083855353   0.011104784   0.005552392   0.003843964 

Para cada una de ellas, observamos cómo se distribuyen el consumo semanal de verdura y el consumo semanal de comida grasa:

tabla = encuesta_salud %>%
  filter(frecuencia_hambre_mensual != 'Dato perdido') %>%
  filter(consumo_semanal_verdura != 'Dato perdido') %>%
  select(consumo_semanal_verdura, frecuencia_hambre_mensual) %>%
  table()
tabla = as.data.frame(tabla)
colnames(tabla) <- c("consumo_semanal_verdura", "frecuencia_hambre_mensual", "q")

order = c("Nunca", "Rara vez", "Algunas veces", "Casi siempre", "Siempre")
type = c("4 o más veces al día", "3 veces al día", "2 veces al día", "1 vez al día", "4 a 6 veces durante los últimos 7 días", "1 a 3 veces durante los últimos 7 días", "No comí verduras ni hortalizas durante los últimos 7 días")


ggplot(tabla) +
 aes(x = frecuencia_hambre_mensual, 
     fill = factor(str_wrap(consumo_semanal_verdura, 20), levels = str_wrap(type, 20)),
     weight = q) +
 geom_bar(position = "fill") +
 scale_fill_hue(direction = 1) +
 coord_flip() +
 theme_minimal() +
 theme(legend.position = "bottom") + 
 theme(legend.title=element_blank()) + 
 labs(x = "Frecuencia de hambre",  y = "Consumo de verdura (proporción)", title = "Consumo de verduras y cantidad de veces con hambre en el mes") + 
 scale_x_discrete(limits = order)

En estos gráficos, se puede observar que aquellos que suelen pasar hambre (siempre o casi siempre) registran un menor consumo de verduras y hortalizas. En los primeros también hay un alto porcentaje de peresonas que comen verduras y/o hortalizas 4 o más veces por día lo que de alguna manera también habla de lo poco balanceadas que son sus dietas.

2) Modelo inicial

La primera alternativa para modelar el peso, es la siguiente:

\[ E(peso) = \beta_{0} + \beta_{1} altura+ \beta_{2} edad+ \beta_{3} genero + \beta_{4} diasActividadF isicaSemanal + \beta_{5} consumoDiarioAlcohol \]

Ajusto el modelo lineal:

modelo_simple = lm(peso ~ altura + edad + genero + dias_actividad_fisica_semanal + consumo_diario_alcohol, data=encuesta_salud)

tidy_modelo_simple <- tidy(modelo_simple, conf.int = TRUE)
tidy_modelo_simple

El valor de \(\beta_{0}\) nos indica que el peso esperado para una mujer de 0 cm de altura, 0 años, que no realiza actividad física y que no consume alcohol es de -68.92 kg.

El coeficiente que multiplica a la altura \(\beta_{1}\) es 0.65. Significa que ante cada cm adicional de altura el peso promedio aumenta en 0,65 kg manteniendo todas las otras variables constantes

El coeficiente que multiplica a generoMasculino es 1.26, de lo que se deduce que el peso esperado para una persona recién nacida de género masculino de 0 cm de altura, que no realiza actividad física ni consume alcohol es de -67,65 (-68.92 + 1.2626).

El valor de \(\beta_{4}\) nos muestra que ante cada día de actividad física adicional el peso promedio disminuye en 0,08kg (manteniendo todas las otras variables constantes).

El valor de \(\beta_{5}\) nos indica que ante cada trago diario adicional de alcohol el peso esperado de una persona aumenta en 0,007kg (manteniendo todas las otras variables constantes).

Sin embargo, no todas las variables son significativas. Tras realizar el test de significatividad individual para los días de actividad física semanal y para el consumo diario de alcohol se observan p-valores mayores a 0.05 por lo que no se rechaza la hipótesis nula (\(\beta_{k}=0\)). También se puede ver que los intervalos de confianza para ambos coeficientes contienen el 0, lo que nos da la pauta de que no son variables útiles para explicar el peso de una persona.

Por otro lado, el resto de las variables (altura, edad y género) sí resultan estadísticamente significativas para explicar el peso de una persona.

glance(modelo_simple)

El \(R^{2}\) del modelo, el porcentaje de variabilidad del fenómeno que el modelo logra explicar, es de 0,35. Considerando el p-valor, se rechaza la hipótesis nula del test de significatividad global y podemos concluir que al menos una de las variables regresoras sirve para explicar el peso de una persona.

3) Modelo categóricas

Luego, se probará con un modelo que incopora el consumo semanal de snacks y una interacción entre el género y la edad, en lugar de actividad física y consumo de alcohol:

\[ E(peso) = \beta_{0} + \beta_{1} altura+ \beta_{2} edad+ \beta_{3} genero + \beta_{4} consumoSemanalSnacks + \beta_{5} genero.edad \]

Se harán algunas modificaciones para que el nivel basal del consumo semanal de snacks sea “no comí comida salada o snacks en los últimos 7 días:

encuesta_salud$consumo_semanal_snacks <- relevel(as.factor(encuesta_salud$consumo_semanal_snacks), ref = "No comí comida salada o snacks en los últimos 7 días")

modelo_categoricas = lm(peso ~ altura + edad + genero + consumo_semanal_snacks + genero * edad, data=encuesta_salud)

tidy_modelo_categoricas <- tidy(modelo_categoricas, conf.int = TRUE)
tidy_modelo_categoricas

En este caso \(\beta_{0}\) nos indica que el peso promedio de una mujer de 0 años, 0 cm, que no consume comida salada ni snacks es de -64.2 kg.

Las distintas categorías de consumoSemanalSnacks (por ejemplo, consumo_semanal_snacks1 a 3 veces durante los últimos 7 días) reflejan la diferencia de peso entre una persona que no consume snacks vs una persona que consume snacks con la frecuencia que refleja la categoría. No todas las categorías de consumo semanal de snacks resultan significativas: aquellas que surgen de consumir snacks todos los días son las que no lo son.

Por otro lado, el coeficiente edad.generoMasculino en este caso refleja que, para los hombres, un año adicional de edad aumenta el peso promedio en 1,22 + 0,39 kg manteniendo constantes el resto de las variables. La variable tiene un efecto significativo para explicar el peso.

Vamos a ver qué porcentaje de la variabilidad explica el modelo:

glance(modelo_categoricas)

Por último, considerando que algunas categorías consumoSemanalSnacks no resultan significativas, se va a realizar un test F para evaluar la significatividad conjunta de la variable para explicar al peso.

tidy(anova(modelo_categoricas))

Considerando el p-valor, se concluye que la variable consumo_semanal_snacks es significativa para explicar el peso de una persona.

A continuación, se va a proponer un nueva definición de las categorías de consumoSemanalSnacks: se van a juntar “consume snacks 1 vez al día”, “consume snacks 2 veces al día”, “consume snacks 3 veces al día” y “consume snacks 4 o más veces al día” en una única categoría “consume snacks todos los días”, y se van a volver a estimar los coeficientes.

encuesta_salud$consumo_semanal_snacks_new <- ifelse(encuesta_salud$consumo_semanal_snacks %in% c("1 vez al día", "2 veces al día", "3 veces al día", "4 o más veces al día"), 
                                                    "Consume snacks todos los días", 
                                                    as.character(encuesta_salud$consumo_semanal_snacks))

encuesta_salud$consumo_semanal_snacks_new <- relevel(as.factor(encuesta_salud$consumo_semanal_snacks_new), ref = "No comí comida salada o snacks en los últimos 7 días")

modelo_categoricas_grouped = lm(peso ~ altura + edad + genero + consumo_semanal_snacks_new + genero * edad, data=encuesta_salud)

tidy_modelo_categoricas_grouped <- tidy(modelo_categoricas_grouped, conf.int = TRUE)
tidy_modelo_categoricas_grouped

NA

Se observa que todas las variables incluídas en el nuevo modelo son significativas (incluídas todas las categorías de consumo semanal de snacks).

A continuación, se va a consultar la variabilidad explicada por el modelo:

glance(modelo_categoricas_grouped)

Con la nueva agrupación no hay una mejora en el \(R^2\) del modelo.

4) Modelos propios y evaluación

A continuación, se sugieren dos posibles modelos adicionales para intentar explicar el peso de una persona:

A: \[ E(peso) = \beta_{0} + \beta_{1} altura+ \beta_{2} edad+ \beta_{3} genero + \beta_{4} consumoSemanalSnacksNew + \\ \beta_{5} consumoSemanalGaseosa + \beta_{6} diasConsumoComidaRapida + \beta_{7} genero.edad + \beta_{8} genero.altura \]

B:

\[ E(peso) = \beta_{0} + \beta_{1} altura+ \beta_{2} edad+ \beta_{3} genero + \beta_{4} consumoSemanalSnacksNew + \beta_{5} frecuenciaHambreMensual + \\ \beta_{6} genero.edad + \beta_{7} genero.altura + \beta_{8} edadConsumoAlcohol . consumoDiarioAlcohol \]

El primer modelo propuesto incorpora una interacción entre género y altura, y el consumo semanal de gaseosa y de comida rápida.

El segundo modelo propuesto, además de la interacción entre género y altura, incorpora la frecuencia de hambre mensual. y también una interacción que consiste en multiplicar la edad en la que comenzaron a consumir alcohol y el consumo diario.

modelo_a = lm(peso ~ altura + edad + genero + consumo_semanal_snacks_new + consumo_semanal_gaseosas + dias_consumo_comida_rapida + genero * edad + genero * altura, data=encuesta_salud)

tidy_modelo_a <- tidy(modelo_a, conf.int = TRUE)
tidy_modelo_a

breve comentario de la significatividad de las variables.

modelo_b = lm(peso ~ altura + edad + genero + consumo_semanal_snacks_new + frecuencia_hambre_mensual + genero * edad + genero * altura + edad_consumo_alcohol * consumo_diario_alcohol, data=encuesta_salud)

tidy_modelo_b <- tidy(modelo_b, conf.int = TRUE)
tidy_modelo_b

breve comentario de la significatividad de las variables.

A continuación, se van a comparar los distintos modelos desarrollados. En primer lugar, se carga el dataset de test en el que vamos a medir la performance de los modelos:

encuesta_salud_test <- read.csv("encuesta_salud_test.csv")

Se crea la variable consumo semanal snacks con la nueva agrupación también para el dataset de test:

encuesta_salud_test$consumo_semanal_snacks_new <- ifelse(encuesta_salud_test$consumo_semanal_snacks %in% c("1 vez al día", "2 veces al día", "3 veces al día", "4 o más veces al día"), 
                                                    "Consume snacks todos los días", 
                                                    as.character(encuesta_salud_test$consumo_semanal_snacks))

Se arma una lista con todos los modelos construídos

models <- list(modelo_simple = modelo_simple, 
               modelo_categorias = modelo_categoricas_grouped, 
               modelo_a = modelo_a, 
               modelo_b = modelo_b)

En el dataset de train, se observa el r2 de todos los modelos

evaluacion_train = map_df(models, glance, .id = "model") %>%
  # ordenamos por R2 ajustado
  arrange(desc(adj.r.squared))

evaluacion_train

breve comentario de los r2 y los r2 ajust.

A continuación, se van a calcular RMSE y MAE tanto para train como para test:

lista_predicciones_training = map(.x = models, .f = augment)
lista_predicciones_testing = map(.x = models, .f = augment, newdata = encuesta_salud_test)

cbind(map_dfr(.x = lista_predicciones_training, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% rename(rmse_train = .estimate) %>% select(modelo, rmse_train) %>% arrange(modelo),
map_dfr(.x = lista_predicciones_training, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% rename(mae_train = .estimate) %>% arrange(modelo) %>% select(mae_train),
map_dfr(.x = lista_predicciones_testing, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% rename(rmse_test = .estimate) %>% arrange(modelo) %>% select(rmse_test),
map_dfr(.x = lista_predicciones_testing, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% rename(mae_test = .estimate) %>% arrange(modelo) %>% select(mae_test)) %>% arrange(rmse_test)

NA

Para definir cual es el mejor modelo, se decide comparar el RMSE en test de los distintos modelo con el objetivo de entender qué tan bien generaliza el modelo en datos nuevos:

Se observa que el mejor modelo es X que tiene el menor RMSE en test. y además, se observa que coincide con que es también el de menor X X X X

5) Diagnóstico del modelo

En la siguiente sección, se observará el cumplimiento de los supuestos del modelo lineal para el modelo inicial.

Los supuestos a analizar son los siguientes:

plot(modelo_simple)

Al utilizar plot() sobre el modelo ajustado, se pueden observar varios gráficos que nos van a permitir analizar los supuestos del modelo lineal:

Residuos vs valores predichos: No parece existir una estructura clara entre los residuos y los valores predichos. Sucede algo similar en el gráfico scale-location.

Normal QQ plot: El extremo superior derecho no se ajusta a la distribución teórica, en este caso la ∼N(0,1), por lo que se deduce que los residuos estandarizados no siguen esa distribución.

Residual vs leverage: Existen algunos puntos con un leverage bastante alto. Vamos a ver cuál es la observación para entender si se trata de un posible outlier:

lista_predicciones_training$modelo_simple %>%
  filter(.hat == max(.hat))

NA

6) Modelo robusto

Por último, se va a leer un nuevo dataset con algunos valores atípicos y vamos a volver a observar la relación entre peso y altura:

encuesta_salud_outliers <- read.csv("encuesta_salud_modelo6.csv")

Observamos el modelo ajustado:

modelo_robusto <- lmRob(peso ~ altura + edad + genero + dias_actividad_fisica_semanal + consumo_diario_alcohol,data = encuesta_salud_outliers)

tidy_modelo_robusto <- tidy(modelo_robusto)
tidy_modelo_robusto

coeficientes estimados:

los coeficientes que multiplican la altura y la edad no cambian sustancialmente, y siguen siendo significativos para explicar el peso de una persona. algo similar sucede con el coeficiente que acompaña al género. Días de actividad física semanal y consumo diario de alcohol siguen siendo no significativas para explicar el peso.

Ahora observamos el porcentaje de la varianza explicada por el modelo:

glance(modelo_robusto)

se observa una baja considerable del \(R^2\) de este nuevo modelo. Pasa de cerca del 0,35 a un 0,28.

A continuación, vamos a observar que tan bien performa el modelo en datos nuevos. Se va a utilizar nuevamente el RMSE y el MAE:

si comparamos el rmse y el mae vs los modelos anteriormente analizados, se puede ver que no se deterioró mucho la predicción del modelo. sí, podría